rm(list = ls())
#working_folder<- "/Users/bgpopescu/Library/CloudStorage/Dropbox/covid_institutions/"
working_folder<-"//Mac/Dropbox-1/covid_paper_replication/"
setwd(working_folder)

#This is to obtain:
#-figure_a13a
#-figure_a13b
#-figure_a13c
#-figure_a13d
#-figure_a13e
#-figure_a13f
#-figure_a13g

library("readxl")
library("reshape")
library("stringi")
library("stringr")
library("dplyr")
library("tidyr")
library("stargazer")
library("corrplot")
library("Hmisc")
library("glue")
library("ggplot2")
library("glue")
library('gridExtra')
library('grid')
library('lfe')
library('broom')
library("broom")
library("ggrepel")
library("lemon")
library("sandwich")
library("multiwayvcov")
library("lmtest")
library("spdep")
library("fastDummies")
library("rlist")
library("vars")
library("conleyreg")
library("tidyverse")
library("geojsonsf")



#################
#Step2: Read CSV#
#################

covid_cases_state <- read_csv("./data/2020_DE_Region_Mobility_Report.csv")
names(covid_cases_state)
covid_cases_state<-subset(covid_cases_state, select = c(date, sub_region_1, 
                                                        retail_and_recreation_percent_change_from_baseline,
                                                        grocery_and_pharmacy_percent_change_from_baseline,
                                                        parks_percent_change_from_baseline,
                                                        transit_stations_percent_change_from_baseline,
                                                        workplaces_percent_change_from_baseline,
                                                        residential_percent_change_from_baseline))
covid_cases_state$date<-as.Date(covid_cases_state$date, "%m/%d/%Y")

covid_cases_state$week <- as.character(strftime(covid_cases_state$date,format="%W"))
covid_cases_state<-subset(covid_cases_state, !is.na(sub_region_1))

covid_cases_state <- covid_cases_state %>% 
  dplyr::group_by(week,sub_region_1) %>% 
  dplyr::summarize(retail_and_recreation_percent_change_from_baseline = mean(retail_and_recreation_percent_change_from_baseline), 
                   grocery_and_pharmacy_percent_change_from_baseline = mean(grocery_and_pharmacy_percent_change_from_baseline), 
                   parks_percent_change_from_baseline = mean(parks_percent_change_from_baseline), 
                   transit_stations_percent_change_from_baseline = mean(transit_stations_percent_change_from_baseline), 
                   workplaces_percent_change_from_baseline = mean(workplaces_percent_change_from_baseline), 
                   residential_percent_change_from_baseline = mean(residential_percent_change_from_baseline))

covid_cases_state$NAME_LATIN<-covid_cases_state$sub_region_1
covid_cases_state<-subset(covid_cases_state, select = -c(sub_region_1))


#####################################
#Summing relevant vars to Land level#
#####################################

data_cov_mob = geojson_sf("./data/data_cov_mob.geojson")
data_cov_mob$pct_manufacturing<-as.numeric(data_cov_mob$pct_manufacturing)

state_weekly<-data_cov_mob%>%
  st_drop_geometry()%>%
  group_by(SN_L, week)%>%
  dplyr::summarize(lag_covid_noneg= sum(lag_covid_noneg, na.rm=T),
                   mean_mobil = mean(sum_mobil, na.rm=T),
                   sum_mobil = sum(sum_mobil, na.rm=T))

state_weekly$post_treat<-ifelse(as.numeric(state_weekly$week)>10, 1,0)


#Step1. Turning some variables to numeric
data_nocovid_nomob<-subset(data_cov_mob, week=="00")

data_nocovid_nomob$gvereine<-as.numeric(data_nocovid_nomob$gvereine)
data_nocovid_nomob$kul<-as.numeric(data_nocovid_nomob$kul)
data_nocovid_nomob$nat<-as.numeric(data_nocovid_nomob$nat)
data_nocovid_nomob$frei<-as.numeric(data_nocovid_nomob$frei)
data_nocovid_nomob$soz<-as.numeric(data_nocovid_nomob$soz)
data_nocovid_nomob$pol<-as.numeric(data_nocovid_nomob$pol)
data_nocovid_nomob$inter<-as.numeric(data_nocovid_nomob$inter)

data_nocovid_nomob$kul<-as.numeric(data_nocovid_nomob$kul)
data_nocovid_nomob$nat<-as.numeric(data_nocovid_nomob$nat)
data_nocovid_nomob$frei<-as.numeric(data_nocovid_nomob$frei)
data_nocovid_nomob$soz<-as.numeric(data_nocovid_nomob$soz)
data_nocovid_nomob$pol<-as.numeric(data_nocovid_nomob$pol)
data_nocovid_nomob$inter<-as.numeric(data_nocovid_nomob$inter)

data_nocovid_nomob$gspo<-as.numeric(data_nocovid_nomob$gspo)
data_nocovid_nomob$gkul<-as.numeric(data_nocovid_nomob$gkul)
data_nocovid_nomob$gnat<-as.numeric(data_nocovid_nomob$gnat)

data_nocovid_nomob$gsoz<-as.numeric(data_nocovid_nomob$gsoz)
data_nocovid_nomob$gpol<-as.numeric(data_nocovid_nomob$gpol)
data_nocovid_nomob$ginter<-as.numeric(data_nocovid_nomob$ginter)
data_nocovid_nomob$pop_60_65<-as.numeric(data_nocovid_nomob$pop_60_65)
data_nocovid_nomob$employment_manufacturing<-as.numeric(data_nocovid_nomob$employment_manufacturing)
data_nocovid_nomob$employment_total<-as.numeric(data_nocovid_nomob$employment_total)

#Summing variables at Land level
state_level<-data_nocovid_nomob%>%
  group_by(SN_L)%>%
  dplyr::summarize(gvereine= sum(gvereine, na.rm=T),
                   gspo= sum(gspo, na.rm=T), 
                   gkul = sum(gkul, na.rm=T), 
                   gnat = sum(gnat, na.rm=T),
                   gsoz = sum(gsoz, na.rm=T),
                   gpol = sum(gpol, na.rm=T),
                   ginter = sum(ginter, na.rm=T),
                   pop_total = sum(pop_total, na.rm=T),
                   pop_60_65 = sum(pop_60_65, na.rm = T),
                   pop_65_75 = sum(pop_65_75, na.rm = T),
                   pop_75_over = sum(pop_75_over, na.rm = T),
                   pct_pop_above_60 = (pop_60_65 + pop_65_75 + pop_75_over)/pop_total*100,
                   log_pop_total = log(pop_total),
                   sum_afd =sum(sum_afd, na.rm=T),
                   sum_voters = sum(sum_voters, na.rm=T),
                   sum_eligible = sum(sum_eligible, na.rm=T),
                   pct_turnout = sum_voters/sum_eligible*100,
                   edu_total = sum(edu_total, na.rm=T),
                   edu_college = sum(edu_college, na.rm=T),
                   pct_college = edu_college/edu_total*100,
                   sum_pop = sum(sum_pop, na.rm=T), 
                   pct_afd = (sum_afd/sum_pop)*100,
                   sum_area = sum(area, na.rm=T),
                   pop_den = (pop_total/sum_area),
                   log_gdp_per_cap = log(sum(gdp_per_cap)),
                   employment_service_sec_narrow = sum(employment_service_sec_narrow, na.rm=T),
                   employment_total = sum(employment_total, na.rm=T),
                   employment_manufacturing = sum(employment_manufacturing, na.rm=T),
                   pct_manufacturing = employment_manufacturing/employment_total*100,
                   pct_service_sec_narrow = employment_service_sec_narrow/employment_total*100)



state_level$bridging_total<-state_level$gspo+ state_level$gkul+ state_level$gnat
state_level$bridging_tot_pop<-state_level$bridging_total/state_level$pop_total*1000
summary(state_level$bridging_tot_pop)
mean_brid<-summary(state_level$bridging_tot_pop)[4]
state_level$bridging_tot_pop_dummy<-ifelse(state_level$bridging_tot_pop>mean_brid, 1, 0)


state_level$vereine_tot_pop<-state_level$gvereine/state_level$pop_total*1000
mean_ver<-summary(state_level$vereine_tot_pop)[4]
state_level$vereine_tot_pop_dummy<-ifelse(state_level$vereine_tot_pop>mean_ver, 1, 0)


state_level$bonding_total<-state_level$gsoz+ state_level$gpol+ state_level$ginter
state_level$bonding_tot_pop<-state_level$bonding_total/state_level$pop_total*1000
mean_bond<-summary(state_level$bonding_tot_pop)[4]
state_level$bonding_tot_pop_dummy<-ifelse(state_level$bonding_tot_pop>mean_bond, 1, 0)

#####################################
#Step1: Reading the state shapefiles#
#####################################


states <-st_read(dsn="./data/shape_files.gdb", layer="state_pt")
states<-st_drop_geometry(states)
states<-subset(states, select = c(NAME_LATIN, CC_1, POINT_X, POINT_Y, NEAR_FID, NEAR_DIST))
states$bfe1<-0
states$bfe1[states$NEAR_FID==4]<-1
states$bfe2<-0
states$bfe2[states$NEAR_FID==5]<-1
states$bfe3<-0
states$bfe3[states$NEAR_FID==6]<-1
states$dist_east_west_brd<-states$NEAR_DIST/1000
states$treat<-NA
states$treat<-0
states$treat[states$NAME_LATIN=="Mecklenburg-Vorpommern" |
               states$NAME_LATIN=="Brandenburg"|
               states$NAME_LATIN=="Berlin"|
               states$NAME_LATIN=="Saxony"|
               states$NAME_LATIN=="Saxony-Anhalt"|
               states$NAME_LATIN=="Thuringia"]<-1


states_poly <- st_read(dsn="./data/shape_files.gdb",
                       layer="gadm36_DEU_1")
states_poly<-subset(states_poly, select = c(NAME_LATIN))
states_poly2<-merge(states_poly, states, by = c("NAME_LATIN"))

##########
#Merging#
#########

states_merged<-left_join(states, state_level, by = c("CC_1"="SN_L"))
state_level<-st_drop_geometry(state_level)
states_poly2_merged<-left_join(states_poly2, state_level, by = c("CC_1" = "SN_L"))
names(states_merged)[names(states_merged)=="CC_1"]<-"SN_L"

final<-left_join(state_weekly, states_merged, by = c("SN_L" = "SN_L"))
final2<-left_join(final, covid_cases_state, by = c("NAME_LATIN" = "NAME_LATIN", "week"="week"))
unique(final2$week)

final2$post_treat<-ifelse(as.numeric(final2$week)>10, 1,0)
final2$post_treat_vereine_tot_pop_dummy<-final2$post_treat*final2$vereine_tot_pop_dummy
final2$post_treat_bridging_tot_pop_dummy<-final2$post_treat*final2$bridging_tot_pop_dummy
final2$post_treat_bonding_tot_pop_dummy<-final2$post_treat*final2$bonding_tot_pop_dummy
final2$lag_covid_pc<-final2$lag_covid_noneg/final2$pop_total*100000
final2$week<-as.numeric(final2$week)+1
final2$week<-as.character(final2$week)
final2$week<-ifelse(as.numeric(final2$week)<10, paste0("0", final2$week), final2$week)

final3 <- dummy_cols(final2, select_columns = 'week')
names_week<-final3 %>% dplyr::select(starts_with("week_"))
list_week_dummies<-names(names_week)


#VEREINE ALL
basic_ver<-c("vereine_tot_pop_dummy", "post_treat", "post_treat_vereine_tot_pop_dummy", "lag_covid_pc")
extended_ver <- c("vereine_tot_pop_dummy", "post_treat", "post_treat_vereine_tot_pop_dummy", "lag_covid_pc",
                  "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                  "pct_college", "pct_pop_above_60", "pct_service_sec_narrow", "pct_manufacturing")

basic_and_dummies_ver<-list.append(basic_ver, list_week_dummies)
extended_and_week_dummies_ver<-list.append(extended_ver, list_week_dummies)



#BRIDGING
basic_brid<-c("bridging_tot_pop_dummy", "post_treat", "post_treat_bridging_tot_pop_dummy", "lag_covid_pc")
extended_brid <- c("bridging_tot_pop_dummy", "post_treat", "post_treat_bridging_tot_pop_dummy", "lag_covid_pc",
                   "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                   "pct_college", "pct_pop_above_60", "pct_service_sec_narrow", "pct_manufacturing")

basic_and_dummies_brid<-list.append(basic_brid, list_week_dummies)
extended_and_week_dummies_brid<-list.append(extended_brid, list_week_dummies)

spec_basic_brid<-as.formula(paste("mean_mobil~",
                                  paste(basic_and_dummies_brid, collapse= "+")))
spec_extended_brid<-as.formula(paste("mean_mobil~",
                                      paste(extended_and_week_dummies_brid, collapse= "+")))

spec_basic_brid_retail<-as.formula(paste("retail_and_recreation_percent_change_from_baseline~",
                                  paste(basic_and_dummies_brid, collapse= "+")))
spec_extended_brid_retail<-as.formula(paste("retail_and_recreation_percent_change_from_baseline~",
                                     paste(extended_and_week_dummies_brid, collapse= "+")))


spec_basic_brid_grocery<-as.formula(paste("grocery_and_pharmacy_percent_change_from_baseline~",
                                         paste(basic_and_dummies_brid, collapse= "+")))
spec_extended_brid_grocery<-as.formula(paste("grocery_and_pharmacy_percent_change_from_baseline~",
                                            paste(extended_and_week_dummies_brid, collapse= "+")))


spec_basic_brid_parks<-as.formula(paste("parks_percent_change_from_baseline~",
                                          paste(basic_and_dummies_brid, collapse= "+")))
spec_extended_brid_parks<-as.formula(paste("parks_percent_change_from_baseline~",
                                             paste(extended_and_week_dummies_brid, collapse= "+")))


spec_basic_brid_transit<-as.formula(paste("transit_stations_percent_change_from_baseline~",
                                        paste(basic_and_dummies_brid, collapse= "+")))
spec_extended_brid_transit<-as.formula(paste("transit_stations_percent_change_from_baseline~",
                                           paste(extended_and_week_dummies_brid, collapse= "+")))


spec_basic_brid_work<-as.formula(paste("workplaces_percent_change_from_baseline~",
                                          paste(basic_and_dummies_brid, collapse= "+")))
spec_extended_brid_work<-as.formula(paste("workplaces_percent_change_from_baseline~",
                                             paste(extended_and_week_dummies_brid, collapse= "+")))


spec_basic_brid_resid<-as.formula(paste("residential_percent_change_from_baseline~",
                                       paste(basic_and_dummies_brid, collapse= "+")))
spec_extended_brid_resid<-as.formula(paste("residential_percent_change_from_baseline~",
                                          paste(extended_and_week_dummies_brid, collapse= "+")))



final2$retail_and_recreation_percent_change_from_baseline
final2$grocery_and_pharmacy_percent_change_from_baseline
final2$parks_percent_change_from_baseline
final2$transit_stations_percent_change_from_baseline
final2$workplaces_percent_change_from_baseline
final2$residential_percent_change_from_baseline


#BONDING
basic_bond<-c("bonding_tot_pop_dummy", "post_treat", "post_treat_bonding_tot_pop_dummy", "lag_covid_pc")
extended_bond <- c("bonding_tot_pop_dummy", "post_treat", "post_treat_bonding_tot_pop_dummy", "lag_covid_pc",
                   "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                   "pct_college", "pct_pop_above_60", "pct_service_sec_narrow", "pct_manufacturing")

basic_and_dummies_bond<-list.append(basic_bond, list_week_dummies)
extended_and_week_dummies_bond<-list.append(extended_bond, list_week_dummies)

spec_basic_bond<-as.formula(paste("mean_mobil~",
                                  paste(basic_and_dummies_bond, collapse= "+")))
spec_basic_int_bond<-as.formula(paste("mean_mobil~",
                                      paste(extended_and_week_dummies_bond, collapse= "+")))




#############################
#Function for calculating SE#
#############################

# Calculate robust confidence intervals
se_robust_cluster <- function(x){
  res<-coeftest(x, vcov = cluster.vcov(x, final3$SN_L, force_posdef=T))[, 2]
  return(res)}

pvalue_robust_cluster <- function(x){
  res<-coeftest(x, vcov = cluster.vcov(x, final3$SN_L, force_posdef=T))[, 4]
  return(res)}

ci_robust_cluster <- function(x){
  res_ci<-coefci(x, vcov = cluster.vcov(x, final3$SN_L, force_posdef=T))
  return(res_ci)}


mo1_bri<-lm(spec_basic_brid, data=final3)
summary(mo1_bri)
mo1_bri_df<-broom::tidy(mo1_bri)
mo1_bri_df$rse<-se_robust_cluster(mo1_bri)
mo1_bri_df$p_rse<-pvalue_robust_cluster(mo1_bri)
mo1_bri_df$std_er_star<-paste("(", round(as.numeric(mo1_bri_df$std.error),3), ")",
                          ifelse(as.numeric(mo1_bri_df$p.value)<0.001,"***",
                                 ifelse(as.numeric(mo1_bri_df$p.value<0.01),"**",
                                        ifelse(as.numeric(mo1_bri_df$p.value)<0.05,"*",
                                               ""))), sep = "")


mo1_bri_retail<-lm(spec_basic_brid_retail, data=final3)
summary(mo1_bri_retail)
mo1_bri_retail_df<-broom::tidy(mo1_bri_retail)
mo1_bri_retail_df$rse<-se_robust_cluster(mo1_bri_retail)
mo1_bri_retail_df$p_rse<-pvalue_robust_cluster(mo1_bri_retail)
mo1_bri_retail_df$std_er_star<-paste("(", round(as.numeric(mo1_bri_retail_df$std.error),3), ")",
                              ifelse(as.numeric(mo1_bri_retail_df$p.value)<0.001,"***",
                                     ifelse(as.numeric(mo1_bri_retail_df$p.value<0.01),"**",
                                            ifelse(as.numeric(mo1_bri_retail_df$p.value)<0.05,"*",
                                                   ""))), sep = "")

mo1_bri_grocery<-lm(spec_basic_brid_grocery, data=final3)
summary(mo1_bri_grocery)
mo1_bri_grocery_df<-broom::tidy(mo1_bri_grocery)
mo1_bri_grocery_df$rse<-se_robust_cluster(mo1_bri_grocery)
mo1_bri_grocery_df$p_rse<-pvalue_robust_cluster(mo1_bri_grocery)
mo1_bri_grocery_df$std_er_star<-paste("(", round(as.numeric(mo1_bri_grocery_df$std.error),3), ")",
                                     ifelse(as.numeric(mo1_bri_grocery_df$p.value)<0.001,"***",
                                            ifelse(as.numeric(mo1_bri_grocery_df$p.value<0.01),"**",
                                                   ifelse(as.numeric(mo1_bri_grocery_df$p.value)<0.05,"*",
                                                          ""))), sep = "")

mo1_bri_parks<-lm(spec_basic_brid_parks, data=final3)
summary(mo1_bri_parks)
mo1_bri_parks_df<-broom::tidy(mo1_bri_parks)
mo1_bri_parks_df$rse<-se_robust_cluster(mo1_bri_parks)
mo1_bri_parks_df$p_rse<-pvalue_robust_cluster(mo1_bri_parks)
mo1_bri_parks_df$std_er_star<-paste("(", round(as.numeric(mo1_bri_parks_df$std.error),3), ")",
                                      ifelse(as.numeric(mo1_bri_parks_df$p.value)<0.001,"***",
                                             ifelse(as.numeric(mo1_bri_parks_df$p.value<0.01),"**",
                                                    ifelse(as.numeric(mo1_bri_parks_df$p.value)<0.05,"*",
                                                           ""))), sep = "")

mo1_bri_transit<-lm(spec_basic_brid_transit, data=final3)
summary(mo1_bri_transit)
mo1_bri_transit_df<-broom::tidy(mo1_bri_transit)
mo1_bri_transit_df$rse<-se_robust_cluster(mo1_bri_transit)
mo1_bri_transit_df$p_rse<-pvalue_robust_cluster(mo1_bri_transit)
mo1_bri_transit_df$std_er_star<-paste("(", round(as.numeric(mo1_bri_transit_df$std.error),3), ")",
                                    ifelse(as.numeric(mo1_bri_transit_df$p.value)<0.001,"***",
                                           ifelse(as.numeric(mo1_bri_transit_df$p.value<0.01),"**",
                                                  ifelse(as.numeric(mo1_bri_transit_df$p.value)<0.05,"*",
                                                         ""))), sep = "")

mo1_bri_work<-lm(spec_basic_brid_work, data=final3)
summary(mo1_bri_work)
mo1_bri_work_df<-broom::tidy(mo1_bri_work)
mo1_bri_work_df$rse<-se_robust_cluster(mo1_bri_work)
mo1_bri_work_df$p_rse<-pvalue_robust_cluster(mo1_bri_work)
mo1_bri_work_df$std_er_star<-paste("(", round(as.numeric(mo1_bri_work_df$std.error),3), ")",
                                      ifelse(as.numeric(mo1_bri_work_df$p.value)<0.001,"***",
                                             ifelse(as.numeric(mo1_bri_work_df$p.value<0.01),"**",
                                                    ifelse(as.numeric(mo1_bri_work_df$p.value)<0.05,"*",
                                                           ""))), sep = "")

mo1_bri_resid<-lm(spec_basic_brid_resid, data=final3)
summary(mo1_bri_resid)
mo1_bri_resid_df<-broom::tidy(mo1_bri_resid)
mo1_bri_resid_df$rse<-se_robust_cluster(mo1_bri_resid)
mo1_bri_resid_df$p_rse<-pvalue_robust_cluster(mo1_bri_resid)
mo1_bri_resid_df$std_er_star<-paste("(", round(as.numeric(mo1_bri_resid_df$std.error),3), ")",
                                   ifelse(as.numeric(mo1_bri_resid_df$p.value)<0.001,"***",
                                          ifelse(as.numeric(mo1_bri_resid_df$p.value<0.01),"**",
                                                 ifelse(as.numeric(mo1_bri_resid_df$p.value)<0.05,"*",
                                                        ""))), sep = "")



#Extended
mo8_bri<-lm(spec_extended_brid, data=final3)
summary(mo8_bri)
mo8_bri_df<-broom::tidy(mo8_bri)
mo8_bri_df$rse<-se_robust_cluster(mo8_bri)
mo8_bri_df$p_rse<-pvalue_robust_cluster(mo8_bri)

mo8_bri_retail<-lm(spec_extended_brid_retail, data=final3)
summary(mo8_bri_retail)
mo8_bri_retail_df<-broom::tidy(mo8_bri_retail)
mo8_bri_retail_df$rse<-se_robust_cluster(mo8_bri_retail)
mo8_bri_retail_df$p_rse<-pvalue_robust_cluster(mo8_bri_retail)

mo8_bri_grocery<-lm(spec_extended_brid_grocery, data=final3)
summary(mo8_bri_grocery)
mo8_bri_grocery_df<-broom::tidy(mo8_bri_grocery)
mo8_bri_grocery_df$rse<-se_robust_cluster(mo8_bri_grocery)
mo8_bri_grocery_df$p_rse<-pvalue_robust_cluster(mo8_bri_grocery)

mo8_bri_parks<-lm(spec_extended_brid_parks, data=final3)
summary(mo8_bri_parks)
mo8_bri_parks_df<-broom::tidy(mo8_bri_parks)
mo8_bri_parks_df$rse<-se_robust_cluster(mo8_bri_parks)
mo8_bri_parks_df$p_rse<-pvalue_robust_cluster(mo8_bri_parks)

mo8_bri_transit<-lm(spec_extended_brid_transit, data=final3)
summary(mo8_bri_transit)
mo8_bri_transit_df<-broom::tidy(mo8_bri_transit)
mo8_bri_transit_df$rse<-se_robust_cluster(mo8_bri_transit)
mo8_bri_transit_df$p_rse<-pvalue_robust_cluster(mo8_bri_transit)

mo8_bri_work<-lm(spec_extended_brid_work, data=final3)
summary(mo8_bri_work)
mo8_bri_work_df<-broom::tidy(mo8_bri_work)
mo8_bri_work_df$rse<-se_robust_cluster(mo8_bri_work)
mo8_bri_work_df$p_rse<-pvalue_robust_cluster(mo8_bri_work)

mo8_bri_resid<-lm(spec_extended_brid_resid, data=final3)
summary(mo8_bri_resid)
mo8_bri_resid_df<-broom::tidy(mo8_bri_resid)
mo8_bri_resid_df$rse<-se_robust_cluster(mo8_bri_resid)
mo8_bri_resid_df$p_rse<-pvalue_robust_cluster(mo8_bri_resid)




######################
#Bridging Time Trends#
######################

df2 <- subset(final3, select = c(mean_mobil,
                               retail_and_recreation_percent_change_from_baseline,
                               grocery_and_pharmacy_percent_change_from_baseline,
                               parks_percent_change_from_baseline,
                               transit_stations_percent_change_from_baseline,
                               workplaces_percent_change_from_baseline,
                               residential_percent_change_from_baseline, week, 
                               bridging_tot_pop_dummy))
df2<-subset(df2, !is.na(bridging_tot_pop_dummy))
df2$type<-ifelse(df2$bridging_tot_pop_dummy==1, "Above Mean", "Below Mean")
unique(df2$type)
df3 <- df2%>%
  group_by(type, week)%>%
  summarise(mean_mobil = mean(mean_mobil, na.rm=T),
            mean_retail = mean(retail_and_recreation_percent_change_from_baseline, na.rm=T),
            mean_grocery = mean(grocery_and_pharmacy_percent_change_from_baseline, na.rm=T),
            mean_parks = mean(parks_percent_change_from_baseline, na.rm=T),
            mean_transit = mean(transit_stations_percent_change_from_baseline, na.rm=T),
            mean_work = mean(workplaces_percent_change_from_baseline, na.rm=T),
            mean_residential = mean(residential_percent_change_from_baseline, na.rm=T))

x<-seq(-10, 42, by = 1)
df3$observation<-NA
df3$observation[df3$type=="Above Mean"] <- x
df3$observation[df3$type=="Below Mean"] <- x

df3$new_type<-NA
df3$new_type[as.numeric(df3$week)<11 & df3$type == "Above Mean"]<-"Above Mean Pre"
df3$new_type[as.numeric(df3$week)<11 & df3$type == "Below Mean"]<-"Below Mean Pre"
df3$new_type[as.numeric(df3$week)>=11 & df3$type == "Above Mean"]<-"Above Mean Post"
df3$new_type[as.numeric(df3$week)>=11 & df3$type == "Below Mean"]<-"Below Mean Post"


df3$observation<-as.character(df3$observation)
selection<-df3$observation[!stringr::str_detect(df3$observation, '-')]

df3$observation[!stringr::str_detect(df3$observation, '-')]<-paste("+", selection, sep = "")
df3$t<-paste("t", df3$observation, sep = "")
df3$t[df3$t=="t+0"]<-"t"


irislabs1<-seq(1,53,6)
irislabs2<-unique(df3[as.numeric(df3$week) %in% irislabs1,]$t)

max_lat<-max(df3$mean_mobil, na.rm=T)
min_lat<-min(df3$mean_mobil, na.rm=T)
max_long<-max(as.numeric(df3$week), na.rm=T)
min_long<-min(as.numeric(df3$week), na.rm=T)
sd_lat<-sd(df3$mean_mobil, na.rm=T)
sd_long<-sd(as.numeric(df3$week), na.rm=T)

coef_int<-round(mo1_bri_df$estimate[mo1_bri_df$term=="post_treat_bridging_tot_pop_dummy"],3)
see_int<-mo1_bri_df$std_er_star[mo1_bri_df$term=="post_treat_bridging_tot_pop_dummy"]
no_int<-paste0(coef_int, "\n", see_int)

bri_mean<-ggplot(df3, aes(x=as.numeric(week), y=mean_mobil, color=type))+
  geom_line()+
  #scale_linetype_manual(name= "", values = c("predicted"="dashed", "actual"="solid")) +
  #guides(linetype = "none")+
  scale_colour_manual(name="", values=c("Above Mean"="red", "Below Mean"="blue"))+
  geom_vline(xintercept = 10, color ="red")+
  #scale_x_continuous(name = "Weeks", breaks=seq(1,55,5))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  labs(x = "Week", y = "Mobility Compared to 2019")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  annotate(geom = "text", min_long + sd_long/5, max_lat-sd_lat/5, label = c(paste("Dif. in Dif.:", 
                                                                                     "\n",
                                                                                   no_int)))
bri_mean<-reposition_legend(bri_mean, 'bottom left')
ggsave(bri_mean, file = "paper/graphs/figure_a13a.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



max_lat<-max(df3$mean_retail, na.rm=T)
min_lat<-min(df3$mean_retail, na.rm=T)
max_long<-max(as.numeric(df3$week), na.rm=T)
min_long<-min(as.numeric(df3$week), na.rm=T)
sd_lat<-sd(df3$mean_retail, na.rm=T)
sd_long<-sd(as.numeric(df3$week), na.rm=T)

coef_int<-round(mo1_bri_retail_df$estimate[mo1_bri_retail_df$term=="post_treat_bridging_tot_pop_dummy"],3)
see_int<-mo1_bri_retail_df$std_er_star[mo1_bri_retail_df$term=="post_treat_bridging_tot_pop_dummy"]
no_int<-paste0(coef_int, "\n", see_int)


bri_retail<-ggplot(df3, aes(x=as.numeric(week), y=mean_retail, color=type))+
  geom_line()+
  #scale_linetype_manual(name= "", values = c("predicted"="dashed", "actual"="solid")) +
  #guides(linetype = "none")+
  scale_colour_manual(name="", values=c("Above Mean"="red", "Below Mean"="blue"))+
  geom_vline(xintercept = 10, color ="red")+
  #scale_x_continuous(name = "Weeks", breaks=seq(1,55,5))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  labs(x = "Week", y = "Mobility Compared to 2019")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  annotate(geom = "text", min_long + sd_long/5, max_lat-sd_lat/5, label = c(paste("Dif. in Dif.:", 
                                                                                  "\n",
                                                                                  no_int)))
bri_retail<-reposition_legend(bri_retail, 'bottom left')
ggsave(bri_retail, file = "paper/graphs/figure_a13b.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)




max_lat<-max(df3$mean_grocery, na.rm=T)
min_lat<-min(df3$mean_grocery, na.rm=T)
max_long<-max(as.numeric(df3$week), na.rm=T)
min_long<-min(as.numeric(df3$week), na.rm=T)
sd_lat<-sd(df3$mean_grocery, na.rm=T)
sd_long<-sd(as.numeric(df3$week), na.rm=T)

coef_int<-round(mo1_bri_grocery_df$estimate[mo1_bri_grocery_df$term=="post_treat_bridging_tot_pop_dummy"],3)
see_int<-mo1_bri_grocery_df$std_er_star[mo1_bri_grocery_df$term=="post_treat_bridging_tot_pop_dummy"]
no_int<-paste0(coef_int, "\n", see_int)

bri_grocery<-ggplot(df3, aes(x=as.numeric(week), y=mean_grocery, color=type))+
  geom_line()+
  #scale_linetype_manual(name= "", values = c("predicted"="dashed", "actual"="solid")) +
  #guides(linetype = "none")+
  scale_colour_manual(name="", values=c("Above Mean"="red", "Below Mean"="blue"))+
  geom_vline(xintercept = 10, color ="red")+
  #scale_x_continuous(name = "Weeks", breaks=seq(1,55,5))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  labs(x = "Week", y = "Mobility Compared to 2019")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  annotate(geom = "text", min_long + sd_long/5, max_lat-sd_lat/5, label = c(paste("Dif. in Dif.:", 
                                                                                  "\n",
                                                                                  no_int)))
bri_grocery<-reposition_legend(bri_grocery, 'bottom left')
ggsave(bri_grocery, file = "paper/graphs/figure_a13c.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)


max_lat<-max(df3$mean_parks, na.rm=T)
min_lat<-min(df3$mean_parks, na.rm=T)
max_long<-max(as.numeric(df3$week), na.rm=T)
min_long<-min(as.numeric(df3$week), na.rm=T)
sd_lat<-sd(df3$mean_parks, na.rm=T)
sd_long<-sd(as.numeric(df3$week), na.rm=T)


coef_int<-round(mo1_bri_parks_df$estimate[mo1_bri_parks_df$term=="post_treat_bridging_tot_pop_dummy"],3)
see_int<-mo1_bri_parks_df$std_er_star[mo1_bri_parks_df$term=="post_treat_bridging_tot_pop_dummy"]
no_int<-paste0(coef_int, "\n", see_int)


bri_parks<-ggplot(df3, aes(x=as.numeric(week), y=mean_parks, color=type))+
  geom_line()+
  #scale_linetype_manual(name= "", values = c("predicted"="dashed", "actual"="solid")) +
  #guides(linetype = "none")+
  scale_colour_manual(name="", values=c("Above Mean"="red", "Below Mean"="blue"))+
  geom_vline(xintercept = 10, color ="red")+
  #scale_x_continuous(name = "Weeks", breaks=seq(1,55,5))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  labs(x = "Week", y = "Mobility Compared to 2019")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  annotate(geom = "text", min_long + sd_long/5, max_lat-sd_lat/5, label = c(paste("Dif. in Dif.:", 
                                                                                  "\n",
                                                                                  no_int)))
bri_parks<-reposition_legend(bri_parks, 'bottom left')
ggsave(bri_parks, file = "paper/graphs/figure_a13d.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)


max_lat<-max(df3$mean_transit, na.rm=T)
min_lat<-min(df3$mean_transit, na.rm=T)
max_long<-max(as.numeric(df3$week), na.rm=T)
min_long<-min(as.numeric(df3$week), na.rm=T)
sd_lat<-sd(df3$mean_transit, na.rm=T)
sd_long<-sd(as.numeric(df3$week), na.rm=T)


coef_int<-round(mo1_bri_transit_df$estimate[mo1_bri_transit_df$term=="post_treat_bridging_tot_pop_dummy"],3)
see_int<-mo1_bri_transit_df$std_er_star[mo1_bri_transit_df$term=="post_treat_bridging_tot_pop_dummy"]
no_int<-paste0(coef_int, "\n", see_int)

bri_transit<-ggplot(df3, aes(x=as.numeric(week), y=mean_transit, color=type))+
  geom_line()+
  #scale_linetype_manual(name= "", values = c("predicted"="dashed", "actual"="solid")) +
  #guides(linetype = "none")+
  scale_colour_manual(name="", values=c("Above Mean"="red", "Below Mean"="blue"))+
  geom_vline(xintercept = 10, color ="red")+
  #scale_x_continuous(name = "Weeks", breaks=seq(1,55,5))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  labs(x = "Week", y = "Mobility Compared to 2019")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  annotate(geom = "text", min_long + sd_long/5, max_lat-sd_lat/5, label = c(paste("Dif. in Dif.:", 
                                                                                  "\n",
                                                                                  no_int)))
bri_transit<-reposition_legend(bri_transit, 'bottom left')
ggsave(bri_transit, file = "paper/graphs/figure_a13e.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



max_lat<-max(df3$mean_work, na.rm=T)
min_lat<-min(df3$mean_work, na.rm=T)
max_long<-max(as.numeric(df3$week), na.rm=T)
min_long<-min(as.numeric(df3$week), na.rm=T)
sd_lat<-sd(df3$mean_work, na.rm=T)
sd_long<-sd(as.numeric(df3$week), na.rm=T)


coef_int<-round(mo1_bri_work_df$estimate[mo1_bri_work_df$term=="post_treat_bridging_tot_pop_dummy"],3)
see_int<-mo1_bri_work_df$std_er_star[mo1_bri_work_df$term=="post_treat_bridging_tot_pop_dummy"]
no_int<-paste0(coef_int, "\n", see_int)


bri_work<-ggplot(df3, aes(x=as.numeric(week), y=mean_work, color=type))+
  geom_line()+
  #scale_linetype_manual(name= "", values = c("predicted"="dashed", "actual"="solid")) +
  #guides(linetype = "none")+
  scale_colour_manual(name="", values=c("Above Mean"="red", "Below Mean"="blue"))+
  geom_vline(xintercept = 10, color ="red")+
  #scale_x_continuous(name = "Weeks", breaks=seq(1,55,5))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  labs(x = "Week", y = "Mobility Compared to 2019")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  annotate(geom = "text", min_long + sd_long/5, max_lat-sd_lat/5, label = c(paste("Dif. in Dif.:", 
                                                                                  "\n",
                                                                                  no_int)))
bri_work<-reposition_legend(bri_work, 'bottom left')
ggsave(bri_work, file = "paper/graphs/figure_a13f.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



max_lat<-max(df3$mean_residential, na.rm=T)
min_lat<-min(df3$mean_residential, na.rm=T)
max_long<-max(as.numeric(df3$week), na.rm=T)
min_long<-min(as.numeric(df3$week), na.rm=T)
sd_lat<-sd(df3$mean_residential, na.rm=T)
sd_long<-sd(as.numeric(df3$week), na.rm=T)



coef_int<-round(mo1_bri_resid_df$estimate[mo1_bri_resid_df$term=="post_treat_bridging_tot_pop_dummy"],3)
see_int<-mo1_bri_resid_df$std_er_star[mo1_bri_resid_df$term=="post_treat_bridging_tot_pop_dummy"]
no_int<-paste0(coef_int, "\n", see_int)


bri_resid<-ggplot(df3, aes(x=as.numeric(week), y=mean_residential, color=type))+
  geom_line()+
  #scale_linetype_manual(name= "", values = c("predicted"="dashed", "actual"="solid")) +
  #guides(linetype = "none")+
  scale_colour_manual(name="", values=c("Above Mean"="red", "Below Mean"="blue"))+
  geom_vline(xintercept = 10, color ="red")+
  #scale_x_continuous(name = "Weeks", breaks=seq(1,55,5))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  labs(x = "Week", y = "Mobility Compared to 2019")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12)) +
  annotate(geom = "text", max_long - sd_long/5, max_lat-sd_lat/5, label = c(paste("Dif. in Dif.:", 
                                                                                  "\n",
                                                                                  no_int)))
bri_resid<-reposition_legend(bri_resid, 'top left')
ggsave(bri_resid, file = "paper/graphs/figure_a13g.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)

